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SUMMARY 


This  report  derives  and  tests  tractable  approximate  formulas  for 
coupling  among  very  low  frequency/low  frequency  (VLF/LF)  waveguide 
modes  that  occurs  at  boundaries  that  separate  regions  having  different 
ground  conductivity.  Although  algebraically  complicated,  these  for¬ 
mulas  are  easily  programmed  and  require  far  less  computer  running  time 
than  the  numerical  mode-coupling  algorithms  used  in  full  wave  computer 
codes . 

To  derive  and  implement  the  formulas,  we  substitute  an  equivalent 
parallel-plate  waveguide  for  the  actual  waveguide  in  the  short  spatial 
interval  that  contains  the  conductivity  boundary.  This  approximation 
is  most  easily  applied  to  daytime  ionospheric  conditions.  We  derive  a 
set  of  coupled  differential  equations  which  have  as  their  first-order 
solution  the  well-known  Wentzel-Kramers-Brillouin  (WKB)  approxima¬ 
tion.  However,  in  contrast  to  the  WKB  solution,  the  method  of  the 
present  report  gives  results  beyond  the  first-order  and  obtains  solu¬ 
tions  to  any  order  necessary  to  achieve  specified  accuracy.  Although 
derived  for  inclusion  in  future  computer  models  of  VLF/LF  worldwide 
atmospheric  noise,  the  mode-coupling  formulas  could  be  used  in  ap¬ 
plications  where  a  large  number  of  propagation  paths  cause  computer 
running  time  to  become  a  problem. 
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PREFACE 


To  assess  the  performance  of  long-wave  communication  links  in 
nuclear  environments,  it  is  necessary  to  calculate  the  effect  of  the 
environment  on  atmospheric  noise  as  well  as  on  the  signal.  Because  no 
data  for  atmospheric  noise  exist  under  nuclear-disturbed  conditions, 
the  Defense  Nuclear  Agency  is  considering  developing  a  computer  model 
that  will  predict  the  noise.  A  major  difficulty  in  implementing  such 
a  model  is  developing  an  accurate,  yet  tractable,  means  of  calculating 
coupling  between  waveguide  modes  over  geologically  complex  regions 
like  Canada  and  Greenland.  Accordingly,  the  present  report  develops 
and  tests  mode-coupling  formulas  suitable  for  inclusion  in  future 
noise  models. 
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SECTION  1 
INTRODUCTION 


Electromagnetic  signals  or  noise  at  frequencies  below  60  kHz 
consist  of  a  number  of  waveguide  modes  that  propagate  in  a  cavity 
bounded  sharply  on  the  bottom  by  the  earth  and  diffusely  on  the  top  by 
the  ionosphere.  The  excitation  and  propagation  of  those  modes  depend 
on  the  electrical  properties  of  the  earth  and  ionosphere.  For  uniform 
paths,  such  as  a  midday  path  over  seawater,  the  modes  propagate  more 
or  less  independently  of  one  another.  However,  on  so-called  mixed 
paths,  along  which  the  properties  of  either  the  earth  or  the  iono¬ 
sphere  change,  the  modes  become  interdependent.  That  process  is 
called  mode  coupling  and  is  strongest  near  transition  zones,  where 
geophysical  changes  are  most  pronounced. 

The  most  common  and  severe  transition  zones  in  the  earth- 
ionosphere  waveguide  are  at  the  boundaries  that  separate  large  regions 
of  dissimilar  ground  conductivity.  Other  transitions  are  at  the 
terminator  and  at  the  edges  of  ionospheric  regions  subjected  to  solar- 
induced  or  nuclear  disturbances.  This  report  mainly  addresses  ground 
conductivity  transitions,  although  the  methods  developed  here  are  also 
applicable  to  other  types. 

Waveguide  transitions  can  be  modeled  as  being  either  gentle  or 
abrupt,  depending  on  whether  they  occur  over  distances  that  are  longer 
or  shorter,  respectively,  than  about  one-sixth  of  the  signal  wave¬ 
length.  Gentle  transitions  cause  only  slight  mode  coupling,  whereas 
abrupt  ones  cause  strong  coupling.  In  the  15  to  60  kHz  band  con¬ 
sidered  here,  one-sixth  of  a  wavelength  never  exceeds  3  km. 

Both  the  WAVEGUID  code  (developed  by  the  Naval  Ocean  Systems 
Center,  San  Diego,  California)  and  the  WAVEPROP  code  (developed  by 
Pacific-Sierra  Research  Corporation)  can  handle  mode  coupling  on  mixed 
paths  having  either  gentle  or  abrupt  transitions.  Data  are  sparse, 
however;  and,  for  lack  of  better  information,  geological  maps  •  sed  for 
this  purpose  show  discontinuous  ground-conductivity  transitio  s  [Wes- 
tinghouse,  1968].  The  codes,  therefore,  treat  conductivity 
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transitions  as  abrupt.  In  fact,  conductivity  transitions  probably 
occur  over  distances  of  several  kilometers  or  more  and  should,  there¬ 
fore,  be  treated  as  gentle. 

Propagation  codes,  which  are  always  expensive  to  run,  are  espe¬ 
cially  costly  for  paths  that  traverse  geophysically  complex  regions, 
such  as  Canada.  Although  that  expense  is  often  manageable,  it  can 
become  overwhelming  if  the  number  of  paths  is  large.  A  case  in  point 
is  prediction  of  low-frequency  atmospheric  noise,  which  is  modeled  by 
calculating  the  energy  that  propagates  from  several  hundred  worldwide 
noise  sources.  The  number  of  propagation  paths  in  that  case  is  enor¬ 
mous,  and  inclusion  of  a  full-fledged,  mode-coupling  matrix  at  each 
conductivity  boundary  is  out  of  the  question.  Simplified  equations 
are  needed  to  make  such  calculations  tractable. 

The  present  report  addresses  the  problem  of  atmospheric  noise 
modeling.  In  order  to  derive  mode-coupling  formulas  that  are  inexpen¬ 
sive  to  compute,  we  make  three  main  assumptions.  First,  because  the 
conductivity  transitions  occur  over  distances  much  shorter  than  the 
earth's  radius,  we  use  a  parallel-plate  waveguide  to  model  the  earth 
and  ionosphere.  That  assumption  allows  calculation  of  mode  parameters 
from  transcendental,  rather  than  differential  equations.  It  is  more 
easily  applied  to  daytime  or  disturbed  conditions,  where  propagation 
is  nearly  isotropic,  than  to  nighttime  conditions,  where  the  iono¬ 
sphere  is  diffuse  and  isotropic  and  therefore  more  difficult  to  model 
as  a  conductivity  half-space.  This  report  does  not  address  nighttime 
propagation.  Second,  we  use  a  generalized  Wentzel-Kramers-Brillouin 
(WKB)  method  to  determine  the  fields.  That  method  is  more  accurate 
than  the  usual  WKB  solutions  (see  for  example  Budden,  1985).  The  main 
validity  criterion  is  that  the  change  in  conductivity  is  gentle  enough 
that  reflections  of  the  modes  can  be  neglected.  The  appendix  shows 
that  criterion  to  be  well  satisfied.  The  third  and  final  assumption 
is  that  relatively  few  modes  need  be  retained.  That  assumption  is 
valid  because  noise  bursts  are  incoherent,  and  phase  need  not  be  kept 
in  mode  summation  (but  is  kept  elsewhere) .  There  is  thus  no  need  to 
keep  many  high-order  modes,  which  are  usually  included  in  computing 
coherent  signals  to  ensure  that  "nulls"  are  calculated  correctly. 


The  new  formulas  connect  long-wave  fields  on  opposite  sides  of 
waveguide  transitions  and  are  in  a  form  that  can  be  used  in  future 
noise  models.  They  require  much  less  computer  time  than  equations' now 
used  in  sophisticated  propagation  codes;  and  they  provide  accuracy  at 
least  commensurate  with — and  probably  better  than — the  accuracy  of  the 
input  data. 
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SECTION  2 

TRANSITIONS  IN  GROUND  CONDUCTIVITY 


In  this  section,  we  discuss  ground  conductivity  maps  used  for 
long-wave  propagation  calculations.  We  also  present  sample  data 
showing  that  substantial  approximations  can  be  made  in  computing  mode 
coupling  at  conductivity  boundaries. 

CONDUCTIVITY  MAPS. 

Long-wave  propagation  does  not  depend  strongly  on  conductivity 
over  regions  where  the  conductivity  exceeds  approximately 
3  X'10“3  s/m,  because  ground  having  such  high  conductivity  reflects 
very  low  frequency/low  frequency  (VLF/LF)  waves  almost  perfectly.  For 
that  reason,  only  slight  mode  coupling  occurs  at  a  boundary  separating 
two  regions  of  high  conductivity.  However,  ground  conductivities 
below  about  lO-^  s/m  strongly  affect  the  excitation  and  propagation  of 
transverse  magnetic  (TM)  signals  in  the  VLF  and  LF  bands,  and  can 
affect  excitation  of  transverse  electric  (TE)  signals.  Such  low 
conductivities  occur  in  Greenland  and  Northeastern  Canada,  as  well  as 
certain  other  regions,  including  parts  of  the  Soviet  Union. 

Figure  1  is  an  example  of  the  conductivity  maps  used  in  long-wave 
propagation  calculations  (Westinghouse ,  1968].  It  shows  VLF  conduc¬ 
tivity  values  throughout  most  of  North  America.  Analogous  maps  exist 
for  virtually  all  regions  of  the  earth.  The  map  divides  Canada  into  a 
number  of  large  regions,  with  uniform  conductivity  assigned  to  each. 
The  regions  are  separated  by  abrupt  boundaries. 

Maps  such  as  that  in  Fig.  1  are  based  on  sparse  data.  There  is 
particularly  little  information  about  remote  areas  like  Northern 
Canada;  for  those  regions,  confidence  in  the  conductivity  values  is 
low  [Westinghouse,  1968].  Although  such  maps  are  useful  for  inferring 
average  ground  conductivity  over  long  propagation  path  segments,  their 
resolution  is  not  good  enough  to  define  the  structure  within  bound¬ 
aries  that  separate  adjacent  regions.  The  boundaries  shown  on  the 
maps  are,  therefore,  not  based  on  data,  and  are  portrayed  as  abrupt 
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simply  for  ease  of  presentation.  Since  other  models  are  unavailable, 
conductivity  maps  like  that  shown  are  used  as  inputs  to  long-wave 
propagation  codes. 

We  assert  that  conductivity  transitions--although  abrupt  on  the 
scale  of  a  continent — are  probably  rather  gentle  on  the  scale  of  a 
wavelength.  Transition  zones  between  conductivity  regions  should, 
therefore,  be  modeled  as  gradual  rather  than  as  sharp  boundaries. 

EFFECTIVE  AND  BULK  GROUND  CONDUCTIVITIES. 

The  electrical  conductivity  of  the  earth  varies  with  depth,  so 
the  effective  ground  conductivity  "seen"  by  a  wave  in  the  earth- 
ionosphere  waveguide  is  an  average  over  about  a  skin-depth  5  from  the 
surface.  Because  the  skin  depth  depends  on  wave  frequency,  the  effec¬ 
tive  ground  conductivity  exhibits  a  weak  frequency  dependence.  The 
map  shown  in  Fig.  1,  for  example,  applies  to  frequencies  from  10  to 
30  kHz . 

Figure  2  shows  skin  depth  versus  frequency  for  several  values  of 
ground  conductivity  a  [Kraichman,  1970].  The  VLF/LF  communication 
bands  are  indicated  by  the  shaded  regions  on  the  right  side  of  the 
figure.  For  normal  ground,  where  a  lies  between  10“^  and  10"^  S/m, 
the  VLF/LF  skin  depths  lie  between  about  30  to  100  m;  for  poorly 
conducting  ground,  where  o  lies  between  about  10-1^  and  10“^  S/m,  the 
skin  depths  lie  between  about  300  to  1000  m. 

Of  all  the  conductivity  boundaries  shown  in  Fig.  1,  those  at 
shorelines  are  expected  to  be  the  most  abrupt.  To  examine  that  be¬ 
havior,  we  use  the  data  given  in  Fig.  3,  which  shows  measured  bulk 
conductivity  versus  distance  Inland  on  the  Olympic  Peninsula  in 
Northwest  Washington  State  [Bostick,  Smith,  and  Boehl,  1977].  Al¬ 
though  not  measured  at  a  low-conductivity  site,  the  data  are  among  the 
few  available  that  illustrate  the  spatial  dependence  of  conductivity 
near  a  boundary. 

Figure  1  shows  that  the  effective  conductivity  of  the  Olympic 
Peninsula  is  10“^  S/m,  which  implies  that  depths  within  several  tens 
of  meters  from  the  surface  contribute  to  the  conductivity.  Figure  3 
shows  in  detail,  on  the  other  hand,  that  at  such  depths  the 
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conductivity  is  about  5  x  lO"^  S/m  at  distances  more  than  20  km  from 
the  shoreline,  and  increases  as  the  shoreline  is  approached.  The 
conductivity  is  presumably  4  S/m  on  the  sea  side  of  the  shoreline,  but 
no  data  were  measured  there. 

We  draw  two  conclusions  from  Fig.  3.  First,  the  inland  effective 
conductivity  on  the  Olympic  Peninsula  (cr  »  5  x  10~3  s/m)  is  within  a 
factor  of  2  of  the  nominal  value  ( a  ~  10~2  s/m)  given  on  the  large- 
scale  conductivity  map.  Considering  the  sparse  data  on  which  the  maps 
are  based,  such  accuracy  is  as  good  as  can  be  expected.  Second,  at 
least  on  the  Olympic  Peninsula,  the  effective  ground  conductivity 
requires  more  than  10  km  to  transition  from  its  shoreline  value  to  its 
inland  value.  Although  those  data  are  site-specific,  we  assert  that 
qualitatively  similar  behavior  occurs  at  other  sites  (where  data  are 
not  available) ,  and  that  even  more  gradual  transitions  probably  occur 
between  the  inland  regions  depicted  in  Fig.  1. 

FULL-WAVE  CALCULATION  OF  FIELDS  AT  TRANSITIONS. 

The  WAVEGUID  [Pappert  and  Shockey,  1972;  Pappert,  Moler,  and 
Shockey,  1970]  and  WAVEPROP  [Field  et  al.,  1976]  codes  calculate 
waveguide  mode  parameters  numerically,  accounting  for  (1)  the  vertical 
inhomogeneity  of  the  ionosphere;  (2)  the  curvature  of  the  earth;  and 
<3)  anisotropy  caused  by  the  geomagnetic  field.  That  calculation 
requires  the  solution  of  nonlinear  differential  equations  for  as  many 
as  10  to  20  modes.  Each  solution  is  subject  to  boundary  conditions  at 
the  earth  and  ionosphere.  That  procedure  uses  substantial  computer 
time ,  even  on  modern  high-speed  machines . 

Despite  their  generality  in  other  respects,  however,  those  codes 
use  equations  that  are  valid  only  when  the  earth-ionosphere  waveguide 
is  uniform  along  the  propagation  path.  Mixed  paths  are  handled  by 
dividing  the  waveguide  into  a  number  of  laterally  uniform  segments, 
calculating  properties  of  the  full  complement  of  modes  in  each  seg¬ 
ment,  and  then  matching  solutions  at  the  segment  boundaries.  For 
example,  analysis  of  a  trans-Canadian  path  across  six  conductivity 
regions  (see  Fig.  1)  would  therefore  require  at  least  six  times  as 
much  computer  time  as  a  path  of  equal  length  over  seawater,  even  if 
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boundaries  between  regions  were  treated  as  abrupt.  For  more  realistic 
conductivity  transitions,  like  that  illustrated  in  Fig.  3,  WAVEGUID 
and  WAVEPROP  must  divide  each  boundary  into  many  uniform  segments--the 
so-called  staircase  fit  to  the  actual  transition.  That  process,  causes 
a  multifold  increase  in  computer  time. 

To  illustrate  the  above  points,  we  use  the  WAVEPROP  code  to 
calculate  fields  near  transitions  of  varying  abruptness.  Figure  4 
shows  an  example  calculated  for  a  transition  from  seawater  (conduc¬ 
tivity  a  =  k  S/m)  to  Greenland  (cr  =  10“^  S/m)  .  The  assumed  frequency 
is  45  kHz,  which  has  a  wavelength  X  of  6.8  km  and  a  reduced  wavelength 
A/2it  of  about  1  km. 

The  solid  line  (Fig.  4)  is  the  electric  field  calculated  for  an 
abrupt  transition;  the  dashed  line  is  the  field  when  conductivity 
changes  linearly  from  4  to  10"^  S/m  over  a  distance  of  1  km.  The 
signal  for  the  abrupt  transition  is  different  from  that  for  the  1  km- 
wide  transition,  so  it  is  necessary  to  account  for  finite  boundary 
widths . 

The  signal  exhibits  peaks  and  nulls  behind  the  transition.  That 
structure  is  caused  by  interference  among  modes  excited  at  the 
shoreline;  we  had  to  retain  20  modes  to  calculate  the  interference 
pattern  accurately.  To  model  the  1  km-wide  transition,  we  had  to 
divide  the  boundary  into  about  30  segments,  most  of  which  were  con¬ 
centrated  between  conductivity  values  of  10"^  to  10“^  S/m.  All  20 
modes  had  to  be  calculated  in  each  of  the  30  boundary  segments ,  so  the 
total  calculation  required  600  (20  x  30)  numerical  solutions  of  the 
field  equations. 

Figure  4  highlights  four  problems  that  must  be  addressed  in 
constructing  a  computer  model  for  atmospheric  noise:  (1)  transitions 
between  geologically  different  regions  must  be  accounted  for  where  the 
conductivity  is  low,  (2)  the  signal  depends  on  the  width  of  the  tran¬ 
sition,  (3)  no  data  exist  that  define  the  correct  transition  widths  in 
most  regions,  and  (4)  full-fledged  calculations  for  gentle  boundaries 
require  an  enormous  number  of  numerical  solutions.  It  is  therefore 
impractical  to  use  WAVEGUID  or  WAVEPROP  directly  in  an  atmospheric 
noise  model  that  accounts  for  propagation  across  regions  of  both  low 
and  laterally  nonuniform  ground  conductivity. 
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NEGLECT  OF  BACKWARD  REFLECTIONS. 

The  well-known  WKB,  or  slowly  varying,  approximation  greatly 
simplifies  propagation  equations  because,  among  other  things,  it 
permits  us  to  neglect  the  signals  reflected  backward  from  gradients  in 
the  propagation  medium.  The  WKB  approximation  can  be  applied  when 
changes  in  the  propagation  medium  occur  over  distances  greater  than 
the  reduced  wavelength  A/2 n  (e.g.,  Budden  [1985]).  If  changes  occur 
over  shorter  distances,  the  transitions  should  be  modeled  as  abrupt, 
and  backward  reflections  should  be  included. 

The  wavelength  A  lies  between  5  and  20  km  in  the  15  to  60  kHz 
communications  bands,  so  the  reduced  wavelength  A  lies  between  0.8  and 
3  km — shorter  than  distances  over  which  large-scale  geophysical  struc¬ 
tures  would  be  expected  to  change,  and  shorter  than  the  lateral  scales 

JL 

indicated,  for  example,  by  Fig.  3.  Those  considerations  lead  to  the 
intuitive  conclusion  that  conductivity  transitions  are  slow  enough 
that  reflections  can  be  ignored,  even  near  shorelines. 

Mode  coupling  calculations  require  application  of  second-order 
formulas,  details  of  which  are  given  in  the  appendix.  We  use  the 
gentleness  of  the  conductivity  transitions  only  to  argue  that  backward 
reflections  can  be  ignored.  The  appendix  also  estimates  the  coeffi¬ 
cient  of  modal  reflection  R  from  the  transition  zone.  For  the  TE 
modes,  R  is  always  small;  for  the  TM  modes,  | R }  <  0.02  in  the  VLF/LF 
range . 


^Although  nongradual  transition  might  occur  in  the  immediate 
vicinity  of  a  shoreline,  where  the  conductivity  changes  from  h  S/m 
(seawater)  to  a  lesser--but  still  high--value  characteristic  of  wet 
ground,  the  important  transitions  to  low  inland  conductivities  occur 
more  slowly,  as  shown  in  Fig.  3. 
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SECTION  3 

MODE-COUPLING  EQUATIONS 


In  this  section,  we  describe  the  model  waveguide  we  use  near 
conductivity  boundaries  and  cite  the  connecting  formulas  based  on  that 
model.  The  appendix  gives  the  derivations  of  those  formulas. 

PARALLEL- PLATE  WAVEGUIDE  MODEL  OF  TRANSITION  ZONE. 

The  major  complexity  in  the  WAVEGUID/WAVEPROP  class  of  codes  is 
caused  by  the  diffuse  ionosphere,  which  requires  that  nonlinear  dif¬ 
ferential  equations  for  the  height  dependence  of  the  wave  admittance 
be  solved  numerically  and  iterated  to  match  boundary  conditions  at  the 
ground.  Those  matching  conditions  give  the  so-called  modal  equation, 
whose  solution  yields  the  attenuation  rates  and  phase  velocities  of 
all  waveguide  modes. 

That  numerical  complexity  can  be  avoided  by  representing  the  true 
ionosphere  by  one  that  is  uniform  and  sharply  bounded.  By  properly 
selecting  the  effective  ionospheric  height  h  and  conductivity  it 
is  possible  to  define  an  equivalent,  sharply  bounded,  earth- ionosphere 
waveguide  that  gives  approximately  the  same  VLF/LF  propagation  as  the 
real  waveguide  [Wait,  1970]. 

The  advantage  of  assuming  a  sharply  bounded  ionosphere  is  that 
the  wave  admittance  is  given  by  a  simple  formula,  rather  than  by  the 
numerical  solution  of  a  complicated  differential  equation.  The  modal 
equation  is  then  a  transcendental  equation,  easily  solved  without 
numerical  integration.  The  disadvantage,  of  course,  is  that  calcula- 
tions  made  using  a  sharply  bounded  model  waveguide  are  less  precise 
than  those  made  using  a  diffuse  ionosphere. 

The  question  is  whether  the  computational  ease  of  the  sharply 
bounded  model  is  worth  the  degraded  precision.  In  most  cases  it  is 
not.  But,  for  the  special  problem  of  computing  mode  coupling  at 
conductivity  transitions,  the  precision  loss  is  minor  and  acceptable. 
The  reason  is  that  the  sharply  bounded  model  need  be  used  only  to 
obtain  connecting  formulas  that  bridge  each  transition,  which  extends 
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over  perhaps  10  to  20  km — an  infinitesimal  distance  on  a  global  scale. 
The  WAVEGUID/WAVEPROP  class  of  codes  can  be  used  on  all  other  portions 
of  the  path. 

Figure  5  is  a  schematic  diagram  of  our  model.  Propagation  is  in 
the  x-direction.  The  ground  conductivity  <j  has  a  constant  value  of 
a  =  «7gj  at  ranges  shorter  than  x  =  0  and  a  different  constant  value 
a  =  <7gH  at  ranges  greater  than  x  =  l .  The  transition  from  a gj  to 
CTgji  occurs  in  the  interval  0  <  x  <  2,  where  a  =  <7g(x)  .  The  abrupt 
boundaries  shown  in  Fig.  1,  and  used  in  most  calculations,  correspond 
to  setting  I  =  0  in  our  model.  The  task  in  the  present  report  is  to 
derive  formulas  that  connect  fields  at  x  <  0  to  fields  at  x  >  2 , 
subject  only  to  the  condition  that  <7g(x)  varies  slowly  enough  to 
permit  us  to  neglect  backward  reflections. 

Outside  the  region  0  <  x  <  £ ,  we  use  full-wave  codes  of  the 
WAVEGUID/WAVEPROP  class,  which  account  for  the  round  earth  and  diffuse 
ionosphere  but  assume  a  uniform  ground.  The  flat  regions  for  0  <  x 
and  x  >  i  indicate  the  region  where  we  match  the  full-wave  modes  to 
the  modes  used  in  the  parallel  plate  waveguide.  In  the  diagram  these 
appear  to  have  a  finite  extent,  but,  in  fact,  we  treat  them  in  the 
calculations  as  if  they  were  infinitesimally  thin. 

CONNECTING  FORMULAS. 

Tractable  equations  that  describe  mode  coupling  at  a  conductivity 
transition  are  derived  in  the  appendix.  Here,  we  present  those  equa¬ 
tions  and  define  certain  parameters.  We  first  write  the  equations  in 
general  notation,  applicable  to  either  polarization,  and  then  special¬ 
ize  the  notation  to  TM  and  TE  modes.  The  equations  given  below  omit 
the  geomagnetic  field  and  are  valid  for  daytime  or  nuclear-disturbed 
ionospheres,  but  not  for  undisturbed  nighttime  ionospheres. 

The  refractive  indices  of  the  ionosphere  and  ground  are: 


14 


Figure  5.  Parallel-plate  waveguide  model  at  conductivity  transition. 
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where  u>  is  the  angular  frequency,  £q  is  the  permittivity  of  free 
space,  and  £r  is  the  relative  permittivity  of  ground.  Note  that  a g, 
£r,  and  hence  ng,  are  functions  of  x  in  the  interval  0  <  x  <  2 . 

The  solutions  of  the  Booker  quartic  [Budden,  1985]  in  the  iono¬ 
sphere  and  ground  are: 


1  +  C 
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and 
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where  Cn  is  the  cosine  of  the  eigenangle  of  the  nth  waveguide  mode  and 
is  found  by  solving  the  modal  equation,  given  below.  To  simplify  the 
notation,  we  occasionally  suppress  the  subscript  n  on  q  and  C;  those 
quantities,  however,  depend  on  the  mode  number. 

Our  goal  is  to  express  each  waveguide  mode  beyond  the  transition 
(x  >  2)  in  terms  of  the  modes  incident  on  the  transition.  For  the  nth 
incident  mode,  we  write: 


-ikSIx 

Fn(x,  z)  =  rnfn(z)  e  n  (X  s  0)  ,  (5) 

where,  k  is  the  wave  number  and,  depending  on  the  polarization,  F  is 
the  y-component  of  the  magnetic  or  electric  fields.  In  region  I  the 
ground  conductivity  is  uniform,  so  the  height-gain  function  f(z)  and 
the  propagation  coefficient  Sj  are  independent  of  x.  Both  are  func¬ 
tions  of  x  in  the  region  0  <  x  £  2,  because  of  the  nonuniform 
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conductivity.  The  factor  Tn  (essentially  the  excitation  factor  of  the 
nth  mode  in  region  I)  is  determined  numerically  using 
WAVEGUID/WAVEPROP  solutions,  which  are  presumed  to  be  available  for 
x  <  0.  The  propagation  coefficient  is  simply 
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backward  reflections, 
the  transition  is: 
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where  An  is  a  normalization  function,  and  is  given  by  Eq.  (19)  for  TM 
and  Eq.  (27)  for  TE  modes.  Equation  (7)  is  the  sought-after  connect¬ 
ing  formula.  We  present  and  discuss  the  terms  of  that  equation  below. 

The  quantity  Qn  is  the  mode-coupling  factor,  and,  as  shown  by  its 
formula,  which  is  given  below,  accounts  for  energy  scattered  into  the 
nth  mode  from  other  modes.  It  also  accounts  for  energy  scattered  from 
the  nth  mode  into  other  modes.  Depending  on  the  particular  situation 
and  the  mode  number,  Q  can  therefore  be  either  positive  or  negative. 
All  terras  other  than  Q  in  Eq.  (7)  involve  only  the  nth  mode,  and 
account  for  changes  in  the  propagation  coefficient  and  height  gain 
caused  by  the  conductivity  contrast  between  and  agl • 

The  mode-coupling  factor  Qn  is  the  crux  of  the  present  report. 

It  is  given  by: 
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The  sum  over  the  index  m  in  Eq .  (9)  accounts  for  scattering  into 
the  nth  mode  from  all  other  modes.  The  term  accounts  for 

scattering  from  the  nth  mode.  In  practice,  it  is  seldom  necessary  to 
carry  modes  higher  than  m  =  5.  The  sum  over  the  index  (i)  in  Eq .  (8) 
accounts  for  various  orders  of  scattering  among  modes.  For  example, 
the  term  in  Eq.  (8)  denotes  single  scattering,  which  involves 

only  energy  transferred  directly  between  the  nth  and  mth  modes.  The 


term  Q^2)  denotes  second-order  scattering,  in  which  energy  is  trans¬ 
ferred  from  all  modes  into  the  mth  mode,  and  then  from  the  mth  to  the 
nth.  The  term  0^3)  denotes  third-order  mode  scattering,  and  so 
forth.  In  practice,  retaining  only  the  single-scattering  term  Q^D  is 
adequate  for  situations  in  which  cTgj  and  OgH  exceed  10"^  S/m.  Higher 
order  scattering  terms  must  be  retained  for  conductivities  lower  than 
10~3  S/m,  although  in  no  case  did  we  find  it  necessary  to  retain  terms 
beyond  (i)  =  4. 

Note  that  g^  in  Eq.  (12)  is  proportional  to  the  term  dfn/dx  lfn 
is  defined  in  Eq.  (20)  for  TM  and  Eq.  (28)  for  TE]  and  is  also  propor¬ 
tional  to  the  lateral  gradient  of  the  ground  conductivity.  For 
uniform  ground,  that  derivative  vanishes  and  there  is  no  mode  cou¬ 
pling,  because  in  that  case  Qn  =  0.  Also,  for  uniform  ground, 

Sn(i)  =  sn(0),  fn(f)  -  fn(0),  and  An(f)  =  An(0)  in  Eq.  (7);  so 
Fjn*(f)  =  F^(0)  ,  as  must  be  the  case  in  the  absence  of  conductivity 
transitions . 

Some  terms  in  Eqs .  (7)  through  (12)  have  different  forms  for  TM 
modes  than  for  TE  modes.  Below,  we  give  those  forms  for  the  two 
polarizations . 

TM  Modes . 

For  TM  modes,  the  amplitude  Fn  is  given  by: 
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ZnH 
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(13) 


where  Zq  =  377  fl  and  Hny  is  the  y-component  of  the  nth  mode  magnetic 
intensity  vector.  The  eigenangle  cosines  Cn  are  the  solutions  of  the 
TM  modal  equation: 


iC^Z^  f  Z^j  cos  kCh  -  +  Z^Z^j  s^n  ^Ch  =  ®  > 


(14) 


where  the  relative  TM  ionospheric  and  ground  impedances  are 
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In  the  transition  region  0  <  x  <  i ,  the  impedance  Zg,  and  hence  Cn,  is 
a  function  of  x.  For  large  conductivities  (a g  >  10“2) ,  where 
| Zg |  «  1,  the  solutions  to  Eq.  (14)  are: 


Cn  *  (n  -  1/2)  ^  ,  n  =  1,  2,  ... 


The  TM  height-gain  function  is 


f  =  cos  kCz  +  — ^  sin  kCz 
n  C 


The  normalizing  function  is 
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And  the  function  f  is  simply 
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In  the  above  equations,  we  have  occasionally  suppressed  the  x- 
dependence  and  the  mode-number  subscript  n. 


TE  Modes. 

For  TE  modes,  the  amplitude  Fn  is  given  by 
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(21) 


where  Eny  is  the  y-component  of  the  nth  mode  electric  vector.  The 
eigenangle  cosines  Cn  are  solutions  of  the  TE  modal  equation: 

ic|zi  +  Zgj  cos  kCh  -  ^1  +  C2  Z^J  sin  kCh  =  0  ,  (22) 

where  the  TE  relative  ionospheric  and  ground  impedances  are 

Z.  =  l/qi  ,  (23) 


and 
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For  large  conductivities,  where  |Zg|  «  1,  the  solution  to 
Eq.  (22)  is: 


C  -  ,  n  =  1,  2 .  (25) 

n  kh 

The  TE  height-gain  function  is 

iq 

f  =  cos  kCz  +  sin  kCz  .  (26) 

n  C 

The  normalizing  function  A^  is 


A  = 


C 


(ikh  +  Z  )  -  q  /C 

O 


(27) 
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And  the  function  C  is 

n 


DETERMINATION  OF  WAVEGUIDE  HEIGHT  h  AND  IONOSPHERIC  CONDUCTIVITY  a t 
FROM  WAVEGUID/WAVEPROP  SOLUTIONS. 

A  major  benefit  of  the  localized  parallel-plate-waveguide  ap¬ 
proach  is  that  it  permits  the  eigenangle  cosines  to  be  found  from 
Eqs .  (14)  and  (22),  the  simple  modal  equations.  The  WAVEGUID/WAVEPROP 
codes,  on  the  other  hand,  must  solve  nonlinear  differential  equations 
and  iteratively  match  boundary  conditions  to  find  Cn.  The  latter 
procedure  requires  much  computer  time  and  may  be  labor-intensive  for 
even  a  skilled  operator.  The  task  remains,  however,  to  find  a  com¬ 
bination  of  effective  waveguide  height  h  and  ionospheric  conductivity 
that  gives  values  of  Cn  that  are  nearly  the  same  as  those  obtained 
using  WAVEGUID/WAVEPROP  in  conjunction  with  realistic  diffuse  iono¬ 
spheres  . 

The  starting  point  for  finding  h  and  is  the  series  of  values 
of  Cn  calculated  with  WAVEGUID/WAVEPROP  in  region  I  (x  <  0— see 
Fig.  5).  By  using  either  the  equation  for  TM  [Eq.  (17)]  or  for  TE 
[Eq.  (25)],  we  can  write: 


h  * 


"  2(Cn+1-Cn) 


Strictly  speaking,  the  value  of  hn  given  by  Eq.  (29)  depends  on  the 
mode  number  n,  but,  in  practice,  h  varies  only  slightly  with  n  (in  our 
experience,  variation  of  h  is  less  than  5  percent.  We  therefore  use 


h  =  i 


where  N  is  typically  a  number  on  the  order  of  5. 


After  finding  h,  we  then  use  the  mode  equations  [Eqs.  (14)  and 
(22)]  to  find  oj_.  Specifically,  we  insert  Cn  (from  WAVEGUID/ 
WAVEPROP) ,  Zg,  and  h  into  those  equations  and  solve  for  Z^,  which 
defines  the  value  of  a^.  As  the  case  for  h,  the  values  of  <7^  calcu¬ 
lated  in  that  manner  depend  only  slightly  on  mode  number  n,  so  we  use 
a  value  averaged  over  the  important  modes . 
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SECTION  4 
TEST  CASES 

In  this  section,  we  present  sample  calculations  that  test  the 
approximate  formulas  of  Sec.  3  against  exact  ones  made  with  the 
WAVEPROP  code.  We  also  discuss  practical  aspects  of  applying  those 
formulas . 


PARALLEL- PLATE  WAVEGUIDE. 

Figure  6  illustrates  how  closely  the  parallel-plate  model 
waveguide  approximates  the  actual  waveguide.  It  shows  ReCn  versus 
ground  conductivity  for  a  frequency  of  45  kHz  and  the  five  lowest 
order  TM  modes.  The  values  of  Cn  (denoted  by  dots  in  the  figure)  were 
calculated  numerically  using  WAVEPROP  and  account  for  earth  curvature, 
the  diffuse  height-profile  of  the  ionosphere,  and  the  geomagnetic 
field.  We  assume  north-south  daytime  propagation  and  a  60  deg  geomag¬ 
netic  dip  angle. 

To  define  the  equivalent  height  and  conductivity  of  the  parallel- 
plate  waveguide,  we  insert  erg  =  10-^  S/m  and  the  corresponding  Cn 
values  from  Fig.  6  into  Eqs.  (14),  (29),  and  (30)  to  find  h  =*  59  km 
and  a i  ~  i0-^  S/m.  Substituting  those  values  of  h  and  a ^  back  into 
Eq.  (14)  gives  the  values  of  ReCn  indicated  by  the  solid  lines  in 
Fig.  6. 

The  agreement  between  the  "exact"  (WAVEPROP)  and  approximate 
(parallel-plate  waveguide)  results  is  excellent  for  all  five  modes  and 
the  entire  range  of  ground  conductivities.*  It  would  be  even  closer 
for  frequencies  lower  than  45  kHz,  which  are  less  affected  by  the 

*The  first  mode  at  10-2  S/m  is  the  only  one  in  error  by  a  substantial 
amount.  This  mode  is  the  most  oblique;  we  expect  it  to  be  affected 
the  most  by  the  curvature  of  the  earth,  which  we  ignore.  That  mode 
is  typically  not  as  important  as  the  second,  since  it  has  higher 
attenuation  and  lower  excitation.  Moreover,  at  high  conductivity 
(<7g  >  10"3  S/m)  there  is  very  little  mode  coupling,  so  the  error  indi¬ 
cated  does  not  matter  for  our  analysis.  However,  this  mode  is  the 
Brewster's  mode  (i.e.,  C  =  Zg)  for  Og  <  10-3  S/m,  where  mode 
coupling  is  important.  In  this  region  the  modes  match  very  well. 
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earth’s  sphericity  and  less  sensitive  to  waveguide  parameters.  Cal¬ 
culations  (not  reproduced  here)  show  equally  close  agreement  for  TE 
modes.  The  modal  equations  [Eqs.  (14)  and  (22)]  may  therefore  be  used 
to  find  Cn(x)  across  the  conductivity  transition.  That  simplification 
reduces  computer  running  time  enormously.  [Note  that  although 
Eqs.  (14)  and  (22)  give  accurate  values  of  |cn|  =  ReCn,  they  do  not 
provide  accurate  values  of  ImCn,  which  govern  the  waveguide  attenua¬ 
tion  rates.  Because  ImCn  «  ReCn,  and  because  the  assumed  transition 
zones  are  too  narrow  {i  «  100  km)  for  attenuation  to  be  important, 
that  inaccuracy  causes  no  problems  for  the  application  addressed  in 
the  present  report.  However,  the  parallel-plate  model  waveguide 
cannot  be  used  to  approximate  long  (>  100  km)  path  segments  where 
attenuation  is  important  and  must  be  calculated  accurately.] 

MODE  COUPLING . 

To  examine  the  dependence  of  mode  coupling  on  the  conductivities 
<7gl  and  ffgH,  we  assume  a  conductivity  transition  where  ag(x)  varies 
according  to  the  formula 


logio  =  '  (31) 

If  we  use  this  logarithmic  variation,  the  results  depend  only  on 
not  on  the  scale  length  H.  In  this  case,  it  is  true  that  mode  cou¬ 
pling  depends  only  on  the  endpoint  values  Cgj  and  ^gii- 

If  nonlogarithmic  transitions  are  used,  the  mode  coupling  depends 
on  transition  shape  as  well  as  on  the  endpoint  values.  That  depen¬ 
dence  is  slight,  however,  and  to  a  first  approximation  it  can  usually 
be  ignored.  Mode  coupling  at  a  boundary  can  therefore  be  calculated 
by  specifying  Ogj  and  &giz  and  using  any  convenient  model  for  Og(x)  in 
the  zone  0  <  x  <  i,  provided  that  the  transition  is  gentle  enough  so 
that  reflections  can  be  ignored.  As  previously  mentioned,  that  near¬ 
independence  of  mode  coupling  on  transition  shape  is  important,  be¬ 
cause  there  are  virtually  no  data  on  how  varies  within  transition 
zones.  fl 

To  illustrate  how  mode  coupling  depends  on  the  conductivity 
contrast  at  a  boundary,  we  assume  a  second-order  (n  =  2)  TM  mode  of 
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unit  amplitude  incident  on  a  transition  zone  where  the  conductivity 
changes  from  a gj  =  1  S/m  to  values  of  dgjj  that  range  from  1  to 
10~5  s/m.  No  other  modes  are  present  for  x  <  0,  so 


t  1 


for  n  =  2  , 
for  n  *  2  , 


(32) 


in  Eq.  (5).  We  further  assume  a  frequency  of  30  kHz. 

Figure  7a  shows  the  amplitude  of  mode  2  versus  &gH,  and 
Figs.  7b-e  show  the  ground-level  amplitudes  of  modes  1,  3,  4,  and  5, 
respectively.  For  ease  of  comparison,  the  scales  used  in  Figs.  7a-e 
are  consistent.  Each  figure  shows  plots  of  several  graphs  correspond¬ 
ing  to  the  various  orders  of  coupling  (values  of  the  index  i)  dis¬ 
cussed  above  in  conjunction  with  Eqs.  (8)  and  (9). 

Figures  7a-e  show  that  when  <7gH  exceeds  about  3  x  10“^  S/m,  the 
amplitude  of  the  second  mode  remains  near  unity  and  the  amplitudes  of 
the  other  modes  are  small.  For  values  of  ^gU  below  3  x  10--*  S/m, 
however,  the  amplitudes  of  modes  1,  3,  4,  and  5  increase  at  the  ex¬ 
pense  of  mode  2.  That  behavior  confirms  our  earlier  assertion  that 
mode  coupling  is  weak  for  high  ground  conductivities,  but  can  be 
strong  for  low  conductivities. 

Another  conclusion  to  be  drawn  from  Figs.  7a-e  is  that  retention 
of  only  first-order  coupling  (i  =  1)  usually  provides  accurate 
results,  although  it  is  necessary  to  retain  as  many  as  three  or  four 
orders  at  the  lowest  values  of  Ogjj.  In  no  case  do  we  find  that 
orders  higher  than  i  =  4  contribute  significantly. 


GENTLE  VERSUS  ABRUPT  BOUNDARIES. 

In  our  final  example,  we  compare  mode  coupling  across  gentle 
boundaries  to  mode  coupling  across  abrupt  boundaries,  such  as  those 
indicated  in  Fig.  1.  We  make  that  comparison  because,  as  discussed  in 
Sec.  2,  mixed-path  calculations  usually  assume  abrupt  boundaries;  and, 
notwithstanding  the  fact  that  the  WKB  approximation  is  invalid  for 
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abrupt  boundaries,  the  WAVEGUID/WAVEPROP  codes  omit  backward-scattered 
modes  for  computational  convenience. 

The  formulas  for  the  gentle  boundary  are  presented  and  discussed 
in  Sec.  3.  For  an  abrupt  boundary,  where  £  =  0,  we  have: 


r 


+ 

n 


1 

2ikS+A+ 
n  n 


(33) 


where  the  superscripted  plus  (+)  and  minus  (-)  denote,  respectively, 
parameters  in  regions  I  and  II.  Backward  reflections  are  also  omitted 
from  Eq.  (33).  This  is  consistent  with  WAVEGUID  and  WAVEPROP.* 

We  again  assume  that  a gj  =  1  and  only  a  second-order  incident 

mode : 


r 

n 


'1  if  n  =  2  , 
0  if  n  *  2  , 


(34) 


but  here  we  use  a  frequency  of  45  kHz,  in  order  to  give  us  a  better 
test  of  the  method.  Figure  8  plots,  as  a  function  of  Ogij,  the  lowest 
five  TM  mode  amplitudes  at  the  beginning  of  region  II,  i.e.,  where 
x  =  0+. 

The  results  for  the  abrupt  and  gentle  cases  are  about  the  same  if 
OgIX  >  10*3  s/m,  again  confirming  that  almost  any  model  will  give  good 
results  if  the  ground  conductivity  is  high  enough.  However,  for 
crgn  <  lO-^  s/m,  the  differences  between  the  two  boundary  types  are 
substantial.  The  first  (n  =  1)  mode  is  the  so-called  Brewster  mode 
which,  although  strongly  excited,  is  very  heavily  attenuated  over 
poorly  conducting  ground.  At  great  distances  beyond  the  boundary,  the 
second  (and  least  attenuated)  mode  is  the  most  important.  Figure  8 


*In  the  sum  of  Eq.  (33),  Zg  - >Zg,  L' Hospital's  rule  gives,  for 

the  m  -  n  term , 


A  S 
n  n 

A+  S+ 
n  n 


r  . 
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bounda 


shows  that,  for  crgU  <  10"^  S/m,  the  abrupt  and  gentle  models  give 
mode-2  amplitudes  that  disagree  by  6  dB .  As  discussed  in  Sec.  2,  we 
believe  the  gentle  boundary  to  be  more  realistic  than  the  abrupt. 
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SECTION  5 
CONCLUSIONS 


.*  c 


In  this  report,  we  have  derived  formulas  for  approximating  VLF/LF 
wave guide -mode  amplitudes  beyond  a  conductivity  boundary  in  terms  of 
the  mode -amplitudes  incident  on  that  boundary.  Although  algebraically 
complicated,  those  formulas  are  easily  programmed  and  require  far  less 
computer  running  time  than  numerical  mode-coupling  algorithms  used  in 
"exact"  computer  codes,  such  as  WAVEGUID/WAVEPROP .  The  formulas  have 
two  desirable  features:  they  are  computationally  simple  and  they 
depend  mainly  on  the  ground  conductivity  values  on  either  side  of  a 
transition,  but  only  slightly  on  the  conductivity  variation  within  the 
transition  itself.  Data  from  available  maps  of  worldwide  ground 
conductivity  can  therefore  be  inserted  directly  into  the  formulas. 

The  mode -coupling  formulas  are  subjected  to  three  areas  of  ap¬ 
proximation,  which  we  believe  are  valid  under  most  circumstances 
encountered  in  practice.  First,  and  most  important,  we  substitute  an 
equivalent  parallel-plate  waveguide  for  the  actual  waveguide  in  the 
short  spatial  interval  that  contains  the  transition  zone.  Second,  we 
ignore  reflections  from  that  zone,  which  requires  that  all  conduc¬ 
tivity  changes  within  the  transition  zone  be  gradual  (occurring  over 
distances  at  least  as  long  as  the  reduced  wavelength  \/2n) .  Third,  we 
neglect  phase  and  use  only  the  magnitudes  of  the  modes  when  performing 
mode  sums.  However,  the  mode-coupling  equations  [Eqs.  (7-12)]  are 
derived  complete  with  phase  terms  that  can,  at  the  discretion  of  the 
analyst,  be  retained. 

The  latter  approximation--neglect  of  phase — is  appropriate  for 
models  of  worldwide  noise  because  such  models  divide  the  earth  into  a 
large  number  of  sections  that  behave  as  equivalent  noise  transmitters. 
Each  noise  transmitter  contains  many  lightning  sources  that  are  uncor¬ 
related  and  noncoherent.  Moreover,  noise  models  involve  average 
rather  than  instantaneous  values,  and  such  averages  smear  out  the 
nulls  that  occur  on  signals  from  coherent  transmitters.  The  main 
effect  of  phase  in  the  mode  sum  is  to  create  nulls;  and  accurate 
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calculation  of  field  strength  in  a  null  requires  retention  of  many 
modes.  Neglect  of  phase  in  the  mode  sum — but  not  elsewhere — 
eliminates  nulls  and  thus  substantially  reduces  the  required  number  of 
modes. 

Although  derived  for  inclusion  in  future  computer  models  of 
VLF/LF  worldwide  atmospheric  noise,  the  mode-coupling  formulas  can  be 
used  in  any  application  where  the  number  of  propagation  paths  is  so 
large  that  the  computer  running  time  becomes  a  problem.  However,  when 
the  formulas  are  used  to  describe  coherent  signals  rather  than  noise, 
the  phase  terms  must  be  retained.  Moreover,  because  retention  of 
phase  necessitates  inclusion  of  more  higher  order  modes,  and  because 
the  effective  height  h  depends  somewhat  on  mode  number  n,  the  mode- 
specific  hn  [Eq.  (29)]  should  be  used  instead  of  the  average  h 
[Eq.  (30)]. 
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APPENDIX 

DERIVATION  OF  EQUATIONS 


INTRODUCTION. 

This  appendix  derives  our  method  for  efficiently  calculating  the 
coupling  of  modes  across  a  region  of  changing  ground  conductivity.  We 
assume  that  the  earth  and  ionosphere  can  be  modeled  locally  as  a 
parallel-plate  waveguide.  That  assumption  leads  to  equations  for  the 
propagation  constant,  which  can  be  quickly  solved  numerically  in  the 
case  of  the  TM  modes  and  analytically  in  the  case  of  TE  modes.  In 
that  way,  we  avoid  the  complicated  mode-coupling  calculations  used  by 
the  more  precise  WAVEPROP  program.  Here,  we  restrict  our  attention  to 
normal  daytime  ionospheric  conditions,  but  the  method  can,  in  prin¬ 
ciple,  be  extended  to  nighttime  conditions. 

We  use  a  generalized  version  of  the  method  developed  by  Wentzel, 
Kramers,  and  Brillouin  for  quantum  mechanical  applications--hence  the 
name  WKB  method.  In  essence,  they  replaced  the  search  for  a  wave- 
equation  solution,  which  is  a  rapidly  varying  function,  by  a  search 
for  a  more  slowly  varying  function.  Under  certain  conditions-- 
typically,  when  the  medium  is  changing  slowly  over  a  wave length- -the 
slowly  varying  function  can  be  assumed  to  be  a  constant.  That  is 
known  as  the  WKB  approximation — or,  stated  more  precisely — the  first- 
order  WKB  approximation.  We  do  not  make  this  assumption,  however,  but 
find  a  coupled  set  of  exact  differential  equations  for  this  more 
slowly  varying  function.  Thus,  our  method  could  be  called  a  general¬ 
ized  WKB  method.  If  we  account  for  both  mode  coupling  and  backward 
mode  reflection  from  the  boundary,  the  equations  are  far  too  complex 
to  be  solved  quickly.  Therefore,  we  make  the  WKB- like  assumption  that 
the  change  in  the  ground  conductivity  is  slow  enough  that  backward 
reflections  can  be  ignored.  At  the  end  of  this  appendix  we  test  that 
assumption  by  making  the  opposite  assumption--that  mode  coupling  can 
be  ignored--and  find  that  the  reflection  is  very  small. 

We  use  Budden's  renormalization  of  the  magnetic  field  (see  Budden 
[1985]).  Thus  the  field  H  used  in  this  appendix  is,  in  fact,  the 
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magnetic  field  multiplied  by  the  impedance  of  free  space.  The  wave 
impedances  are  therefore  dimensionless  and  should  be  multiplied  by 
377  n  to  put  them  into  MKS  units.  The  above  also  removes  the  distinc¬ 
tion  between  the  magnetic  and  electric  fields.  In  this  appendix,  we 
often  combine  the  standard  TE  (Ey,  Hx,  Hz)  or  TM  field  components  (Hy, 
Ex,  Ez)  into  a  new  field  called  F.  That  notation  allows  us  to  derive 
one  set  of  equations  that  is  applicable  to  both  TE  and  TM  modes. 

We  also  assume  that  the  impedance  of  the  earth  is  independent  of 
the  mode  parameter  (Cn) .  As  shown  later,  that  approximation  is  tan¬ 
tamount  to  assuming  qg  is  independent  of  Cn.  The  definition  of  qg  is 
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The  mode  dependence  is  second  order  in  Cn/ng.  If  we  make  the  restric¬ 
tion  that  ( Cn |  s  0.5,  then  the  largest  error  is  at  very  low  conduc¬ 
tivity  (erg  ~10-5  S/m),  where  it  is  about  2  percent.  That  approxima¬ 
tion  is  not  critical  to  the  analysis,  but  allows  analytic  forms  to  be 
derived  for  the  TE  case.  Without  that  approximation,  integrals  over 
height  from  0  to  «,  would  have  to  be  changed  to  integrals  from  -<*>  to 
®,  and  some  other  analytic  forms  will  change. 

Many  quantities  used  here  are  functions  of  distance  from  the 
transmitter  (denoted  by  x) ,  height  above  the  ground  (z) ,  and  mode 
number  (n) .  When  it  does  not  cause  confusion,  we  omit  the  function 
parameters  or  mode  number  subscript  in  order  to  avoid  visual  com¬ 
plexity. 


MODEL. 

We  treat  the  earth  as  a  conducting  half-space  below  z  =  0,  and 
the  ionosphere  as  a  conducting  half-space  above  z  =  h.  The  index  of 
refraction  for  the  earth  and  the  ionosphere  is,  respectively, 


2 

n 

g 


a 

i  -J5- 
W£0 


(36a) 


40 


* 


(36b) 


Here,  we  assume  that  a g  >  10"^  S/m,  but  we  make  no  assumption  (momen¬ 
tarily)  about  <7£.  We  also  assume  that  the  conductivity  of  the  earth 
is  constant  for  x  <  0  and  x  >  A ,  but  that  it  changes  in  the  region 
0  <  x  <  i .  For  our  purpose,  the  conductivity  is  considered  to  be 
continuous.  We  assume  that  the  boundary  is  far  enough  from  the  trans¬ 
mitter  that  the  total  field  is  made  up  of  a  sum  of  fields  of  in¬ 
dividual  waveguide  modes.  For  example,  the  electric  field  is 


E  =  >  E 

y  Z-rf  yn 

n 


MAXWELL'S  EQUATIONS. 

By  using  Budden's  renormalization  we  can  show  that  Maxwell's 
equations  for  the  TE  case  become: 


.  dE 

H  =  — - ^ 

x  ik  dz  ’ 


.  3E 

H  =  -  — - £ 

z  ik  3x  ' 


ikn2(x)  -  E  +  -j-  H  =  0 

lk  .  2  y  3x  z 
3  z 


For  the  TM  case,  Maxwell's  equations  become: 


ikn  E  =  -  -r-2  , 
x  3z 
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(42) 


(43) 


FIELDS  WHEN  ag  CHANGES. 

To  illustrate  our  technique,  we  use  Fy,  the  y  component  of  the 
field,  in  the  region  where  ground  conductivity  is  not  changing,  i.e., 
x  <  0.  (Note  that  Fy  =  Ey  in  the  TE  case  and  Fy  =  Hy  in  the  TM.)  The 
form  of  Fy  is  well  known;  i.e.f 

_  \  '  -  .  .  -ikSx  ,  _  ikSxl  .... 

Fy  ’  L,  £n(Z)  lrn  '  +  Rn  6  J  '  <44) 

n 

Here  fn(z)  is  the  height  gain  function,  Tn  e“*-kSx  represents  the 
forward-moving  wave,  and  Rn  e^Sx  represents  the  backward-moving  wave, 
i.e.,  that  part  of  the  wave  reflected  from  the  boundary.  We  normalize 
the  height  gain  to  one  on  the  ground,  so  that  the  excitation  (Tn)  and 
reflection  (Rn)  factors  represent  the  forward  and  backward  field  on 
the  ground  at  x  =  0. 

In  the  region  where  the  ground  conductivity  is  changing  we  need  a 
form  that  reduces  to  Eq.  (44).  Thus,  for  x  >  0, 


Fy  *  £  £n<Z’  X)  [An(x)  +  Vx)]  '  <45) 

n-1 


Here,  An(x)  represents  the  forward-moving  wave,  Bn(x)  represents  the 
backward-moving  wave,  and  fn(z,  x)  is  the  height  gain,  which  is 
dependent  on  x.  Here,  An  and  Bn  are  arbitrary  functions  of  x.  We 
place  restrictions  on  the  height  gain  term  fn,  and  determine  what  form 
An(x)  and  Bn(x)  must  have  to  satisfy  Maxwell's  equations. 


42 


The  restriction  on  fn  is  that  locally  it  must  satisfy  the  wave 
equation 

2 

■—  +  k2  q2 (x ,  z)  f  (z,  x)  =  0  ,  (46) 

.  dz  J 

where  q2  =  n^(x,  z)  -  1  +  C^. 

We  discuss  the  height-gain  term  fn  in  more  detail  below;  there  we 
show  that  fn  is  orthogonal  in  the  sense  that 

f  (z,  x)  f  (z,  x)  dz  =  A  (x)  6  .  (47) 

n  m  n  nm 

The  boundary  conditions  satisfied  by  fn  determine  the  modal  equation 
that  gives  the  value  of  Cn.  Thus,  the  information  about  the  ground 
conductivity  is  contained  in  fn,  so  we  can  restrict  ourselves  to  the 
region  within  the  waveguide  (0  <  z  <  h) . 

Equations  (39)  and  (42)  can  be  written: 


F 

z 


ii, 

ik  3x  ‘y  * 


(48) 


(So  Fz  =  Hz  in  the  TE  case  and  Fz  =  -Ez  in  the  TM  case.)  The  form  of 
Fz  that  is  consistent  with  Eq .  (45)  and  that  reduces  to  the  correct 
form  for  x  <  0  is 


We  combine  E.-;s.  (40)  and  (43)  (in  the  waveguide)  to  give: 


_6_ 

<3x 


F  . 
z 


(49) 


(50) 
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A  » 


IkS  B 
n  n 


K  (A«  *  +  5  E 


B  1  +  r—  K  =0 
m  S  mn 
n 


Now  we  define  A  and  B  in  terms  of  a  and  b,  which  are  functions  that 
vary  more  slowly  with  respect  to  x,  and  we  let 


a  (x)  -lk  J\_  S  (x  )  dx 
.  n  J0  n 

A  =  — — — —  e 

n  ys  (x) 


b  (x)  ik  S  (x' )  dx' 
n  J0  n 


B  = 
n 


ysn<x) 


Thus  Eqs.  (54)  and  (55)  become 


i  S'  2ik  r*  S  (x')  dx' 

,  In,  J  0  n 

a'  -  x  b  e 
n  i.  S  n 
n 


S  +  S 
n  ra 


i  ys~^~ 

1  I  I  n  m 


-lk  iS  [Sm(x,)-Sn(x'>]  dx' 


S  -  S  lk  J*  S  (x')  ♦  S  <x')|  dx' 

+  -S - -  be  0  l*  n  J  ■  K  -  0 

rz — —  m  mn 
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b' 

n 


!s;  -2lk  ^  Sn(x-)  dx’ 

2  S  an  e 
n 


5  E 


m~l 


S  -  S 
n  m 


J  S  S 

n  m 


a  e 
m 


■‘kJ?  [Sm(x')+Sn<x,)]  dx' 


s  +  s 

n  m 


J  S  S 

n  m 


b  e 
m 


lk/x  [sB(x')-Sn(*->]  dx’l 


K  =  0  . 
mn 


(59) 


Equations  (56)  and  (57)  would  be  the  WKB  approximation  if  an(x)  and 
bn(x)  were  constants.  Thus,  Eqs.  (58)  and  (59)  become  the  second- 
order  WKB  approximation.  Up  to  this  point,  we  have  made  no  approxi¬ 
mations.  On  pp.  66-68  we  show  by  direct  calculation  that  the  reflec¬ 
tion  term  is  very  small.  Therefore,  we  can  neglect  reflections.  We 
can  now  set  bm  =  0  in  Eq.  (58). 

Removing  the  nth  term  from  the  summation  in  Eq.  (58)  yields: 


a'  +  a  K  + 
n  n  nn 


CO 

E 

S  +  S 
n  m 

Js  S 

m=l 

nvn 

l  m  nj 

a  e 
m 


-ik/J  [smCx')-Sn<x'>]  dx'  k 

mn 


=  0 


(60) 


That  reduces  to 


a  very  simple  form  if  we  let 


a  (x) 
n 


“n<X) 


Recall  that,  from  Eq.  (53), 


(61) 
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K  (x)  -  y- 
nn  A 


n 


f'  (z)  f  (z)  dz 
n  n 


However,  from  Eq.  (47)  we  see  that 


a;(x) 


f'  (z)  f  (z)  dz 
n  n 


So , 


K  =  7T  > 
nn  2A  n 
n 


and  hence , 


K  dx 
nn 


-  5  [An<;O/An(0)] 


Thus  Eq.  (61)  is 


Vr  (0) 

FTix)  • 


So,  from  Eq.  (60)  it  follows: 


„;(x)  - 


E 


g  (x)  «(x)  > 
°nm  m 


m=»l 

mHn 


where  for,  m  *  n, 
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(62) 


(63) 


(64) 


(65) 


(66) 


(67) 


g  (X) 
run 


1 

2 


S  +  S 
n  m 


Js S 

m  n 


A(0)  A  (x) 
m  n 


A  (x)  A  (0) 
m  n 


-ik  f*  fs  -S  1  dx’ 

X  e  0  <■  “  "J  K 


mn 


We  show  below  that  the  function  Kmn  (when  m  *  n)  has  the  form: 


k  44 


1 


mn  A  ik  _2  „2  <3x 

n  C  -  C 
n  m 


TZ  fOO  . 


where  f(x)  =  qg(x)  for  the  TE  case,  or  f(x)  =  Zg(x)  for  the  TM. 
Thus  putting  Eq.  (69)  into  Eq.  (68)  we  get 


g  (x)  = 

nm 


-ik  f*fs  -S  l dx' 
.  J0I  m  n 

1  e  v  ' 


2ik  JS  r  fS  -  S  ] 
n  mi  m  nj 


A  (0) 
m 


r  (*) 


A  (0)  J A  (x)  A  (x) 
n  n  m 


Here,  we  used 


C2  -  C2  -  1  -  C2  -  [l  -  C21  =  S2  -  s2  =  [s  -  S  Ifs  +  S  ] 
nm  nj  mnimnjimnj 


We  can  integrate  Eq.  (67)  directly,  then: 


a  (x)  =  a  (0)  + 
n  n 


E 


rx 


m»l 

m*n 


g  (x')  a  (x' )  dx' 
nm  m 
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(68) 


(69) 


(70) 


(71) 


(72) 


*- 

v 

» t 


o, 

t 

5 

K 


l 


i*  ***  **  ** 
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or,  in  general: 


G(i)(x) 

nm 


E  •  E 


/•X 


dx’  Sn 


dx"'  gkn,<*"’> 


(78) 


E 


eX 


d*'  S„ !<*'>  ^  <*'> 


(79) 


For  x  =  0,  and  ignoring  reflections,  we  have,  from  Eq.  (4-4) 


f  -  V*  r  f  (z)  . 

y  n  n 

n-1 


(80) 


At  x  =  0,  we  have,  from  Eq.  (45): 


E  £n(z’  °>  V0) 


n=l 


(81) 


These  must  be  equal.  So,  since  from  Eq .  (56)  and  Eq.  (61) 


A  (0)  =  a  (0)/JSJ0)  =  a  (0)/iSJ0)  , 
n  n  n  n  u 


(82) 


we  see  that 


Qn<°>  -  rn 


(83) 


That  is  the  initial  condition  on  a. 
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HEIGHT  GAIN  AND  MODAL  EQUATIONS. 

TE  Case. 

In  the  TE  case,  we  write  the  height-gain  function  as  en(z) . 
Recall  that  it  is  a  solution  to  Eq .  (46).  The  boundary  conditions 
require  that  en  and 
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(88) 


e  ■  7T  r  e 
n  lk  dz  n 

be  continuous  across  the  boundaries  at  z  =  0  and  z  =  h.  Thus,  we  can 
define  an  impedance  function: 

e  (z) 

Z(z)  -  — -  .  (89) 

Therefore,  the  boundary  conditions  are 

2(0)  •  -t  .  Z  ,  (90) 

qg  g 

and 

Z(h>  -  -  -f-  -  -Z.  .  (91) 

"l  1  1 


The  solutions  of  Eq.  (46)  that  satisfy  Eq.  (90),  Eq.  (91),  and  the 
radiation  condition  (i.e.,  that  the  field  must  go  to  zero  at 
infinity),  and  have  en(0)  =  1  are: 


e 

n 


\  K  -  L1  -ickh  -^ik<z-h) 

Z  CZ.  +  1  e  e 

S  i 


cos  Ckz  +  i  sin  Ckz 


iq  kz 
e  g 


for  z  >  h  , 


for  0  <  z  <  h  , 


for  z  <  0  . 


(92) 


(93) 


(94) 


Recall  that  the  solution  to  Eq.  (46)  is  a  superposition  of  up-  and 
down-going  waves  whose  amplitudes  are  constants  to  be  determined  by 
the  boundary  conditions.  Thus,  in  the  three  regions  there  are  six 


52 


constants.  Two  of  them  are  determined  by  the  radiation  condition  and 
one  by  the  normalization  condition.  In  reality,  the  boundary  condi¬ 
tions  at  z  =  0  and  z  =  h  determine  four  more  constants,  but  there  are 
only  three  to  be  determined.  That  fact  places  a  restriction  on  the 
value  of  C,  which  is  expressed  by  the  modal  equation: 


CZ  +  1  CZ.  -  1 

g  _  r  -2iCkh 

CZ  -  1  '  CZ.  +  1  e 
S  i 


Note  that  we  can  write  Eq.  (95)  as 


:K +  zJ 


cos  Ckh  —  1  +  C  Z.Z  sin  Ckh 


The  last  form  is  particularly  easy  to  solve  numerically.  When  a g  is 
very  large,  Zg  is  very  small,  so: 


iCZ.  cos  Ckh  -sin  Ckh  . 
l 


To  zeroth  order  in  Z^,  nw/kh,  while  to  first  order 


C(1)=  C(0) 


n  ~  n  ll  -  iZ^/kh  I 


We  solve  Eq .  (96)  numerically  using  the  Newton-Raphsori  method  starting 
with  C  =  C^1) . 

TM  Case. 

In  this  case,  the  height-gain  function  is  written  as  hn(z)  and 
the  boundary  conditions  require  that  hn  and 


■ _ 1_  J_  , 

\  "  . .  2  dz  n 

ikn 
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be  continuous  across  the  boundaries.  Thus, 


Z  (z) 
n 


h  (z) 
n 

h  (z) 
n 


must  also  be  continuous.  The  boundary  conditions  are  then: 


Z(0) 


and 


Z(h) 


Z 


i  ' 


This  gives: 


C  -  Z 

cTz. 

i 


g  -iCkh 
b  e 


-iqik(z-h) 

e 


for  z  >  h  , 


hn(z,  x)  = 


cos  Ckz  +  i  sin  Ckz 
0 


iq  kz 

,  S 


for  0  <  z  <  h, 


for  z  <  0  , 


and  the  modal  equation 


C  +  Z  C  -  Z 
_ _ _ _ g  _  _ i  -ZiCkh 

c-z  ~  c  +  z  e 

g  i- 


Equation  (106)  may  also  be  written: 


(100) 

(101) 

(102) 

(103) 

(104) 

(105) 

(106) 
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icj^Z^  +  zj  cos  Ckh 


(c2  +  hzs) 


sin  Ckh 


(107) 


The  zeroth  order  solution  for  small  Zg  is  then 


r<°)  _  L  I)  «■ 

n  (  2  J  kh  ’ 


(108) 


and  the  first  order  is 


C(1)  =  C(0)  fl  -  - - - 

n  n  1  +  iZ.kh 

L 


(109) 


ORTHOGONALITY  OF  MODES. 

To  prove  the  orthogonality  of  fn(z,  x)  we  note  that,  since  fn  is 
a  solution  of  Eq.  (46),  we  can  write: 


f  €  e  a  ,  x  ,  2 

r,  — t  f  -  f  — -  f  +  k 
m  _  2  n  n  „  2  m 
a  z  dz 


fc2  -  C2|  f 

1  n  mj  n 


f  =  0  . 
m 


(110) 


If  n  *  rn,  then 


f  f  dx  = 

ii  m 


!fc?  -  <1 

V  11  "'J 


~ff  f  "f  4-  f  1  dz 

oz  (  m  3z  n  n  dz  m  J 
n 


(1U) 


The  radiation  condition  ensures  that  fm  ->  0  as  z  ->  m.  The  fact  that 
fm  and  5fn/9z  are  continuous  across  the  boundary  at  z  =  h  means  that 
there  is  no  contribution  to  the  integral  at  z  =  h.  If  the  lower  limit 
of  the  integral  goes  to  -«>,  then  clearly,  Eq .  (Ill)  would  reduce  to 
zero.  However,  using  the  current  limits,  since  fn(0)  -  1,  Eq .  (Ill) 
reduces  to: 
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jyi/lfu  sn  l  /v  if-  try  iru 


~  f'  +  2k2q'f  +  k2q2f'  -  0  . 
„  2  n  n  n  n  n 

dz 


(117) 


MulCiplying  Eq .  (46)  for  mode  m  by  f and  subtracting  Eq .  (117)  multi¬ 
plied  by  fm  yields: 


m  dz2 


2 

f'  -  f'  f  +  2k2q'q  f  £  +  k2fq2  -  Q?)  f'f  =  0  .  (118) 

n  n  2  m  n  m  J/n  J  n  m 


Assuming  n  r*  m: 


f- 

n 

J0 


f  dz  = 
m 


2  21 

C  -  C 
n  in  I 


dz  j-  [f  f-  £'  -  £'  f.  I 
dz  1  m  dz  n  n  dz  ml 


(119) 


which  gives  (after  due  attention  to  the  boundary  conditions): 


K  -  Cn) 


af' 

af 

a 

f ,  m 

dz 

n  dz 

* 

. 

'  z=0 

(120) 


If  we  write 


f  (z,  x)  =  cos  Ckz 


sin  Ckz  , 


(121) 


then  Eq.  (120)  becomes 


K  -  +  T 
nm  A 


run  a  r  0  ^ 

ffl  ikfc2  -  C2] 

l  in  nj 

where  £(x)  =  qg  in  the  TE  case  and  f(x)  =  Zg  in  the  TM. 


(122) 
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MODE  COUPLING  AT  SHARP  BOUNDARY. 

Using  an  adaptation  of  the  method  of  Pappert  and  Shockey  [1972], 
we  can  compare  our  method  with  the  results  from  methods  using  a  sharp 
boundary . 

First  we  write  Fyn  and  Fzn  as  the  y  and  z  components  of  the  field 
of  the  nth  mode.  In  the  TE  case,  we  take  Fz  =  H2  and  in  the  TM  case, 
Fz  =  -Ez .  If  we  write: 


L  = 


1+2  2 

k  3z 


1  ’ 

0 


and 


u 

n 


(123) 


(124) 


then  Maxwell's  equations  can  be  written: 


_L  _3_  - 

ik  5x  Un 


U25) 


If  n|  is  not  a  function  of  x,  then  x  appears  only  in  an  exponential 
teim,  as  for  example; 


-ikS  x 

u  =*rf(z)e  n  (126) 

yn  n  n 


So  Eq .  (125)  becomes: 


A  — V  “+ 

Lu  S  u 
n  n  n 


(127) 


Now  we  define  vv  to  be  eigenvectors  of  ,  the  transpose  of  L,  and 
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introduce  the  following  bracket  notation,  for  arbitrary  vector  func¬ 
tions  r(z)  and  s(z) 


<r,  s>  * 


r(z)  •  s(z)  dz  . 


(128) 


Now  it  is  clear  that 


<v  ,  Lu  >  =  <L  v  ,  u  > 
m  n  m  n 


If  we  let  the  eigenvalue  of  vm  be  so 


'  f  -»  -» 

L  v  =  v  v  , 
m  mm 


then  Eq.  (129)  becomes: 


S  <v  ,  u  >  =  v  <v  ,  u  > 
n  in  n  m  in  n 


Finally  we  see: 


It  is  clear  that 


v  = 
m 


F_1 


ZlT! 


F 

l  ym 


JO  , 
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(129) 


(130) 


(131) 


D  —  1/  ciiiu  ^  —  V  li  LI  r  Hi 

n  ii  m  n 


(132) 


(133) 


*i.v 


<v  .  u  > 
m  11 


F 

pOO 

(F  ,  F  "j  • 

yn 

dz  = 

fF  F  +  F  F  1 

(  zm  ymj 

F 

i  ym  zn  zm  ynj 

J0 

k  zn 

J 

0 

(134) 


Therefore,  since  Fzm  =  snfyn  [Eq.  (127)], 


-4-  M- 

<v  ,  u  > 
m  n 


fS  +  S  ) 
{  n  mj 


F  F  dz 

ym  yn 


(135) 


2  S  A  6 
m  m  nrn 


(136) 


For  the  case  of  waves  on  both  sides  of  the  boundary,  let  represent 
the  wave  before  the  boundary  and  represent  the  wave  after  it.  Mow 
the  boundary  condition  (that  Fy  and  Fz  be  continuous)  is  expressed  as 


E  %  ’  E 


-►+ 

u 

n 


(137) 


n=l 


n=l 


If  we  multiply  by  vj  and  integrate  from  0  to  we  find  that: 


E 

n=l 


-*+ 

<v  ,  u  >  = 

m  n 


E 

n=l 


'4"f  “►+ 

<v  ,  u  > 
rn  n 


(138) 


From  Eq.  (134)  it  is  clear  that: 


<U  ,  V  > 

m  n 


■  lS*  +  Sn) 


F+  F‘  dz 
ym  yn 


(139) 


s+  +  S~  X 
hi  nj  mn 


(140) 
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yu  .  \T\  i  -Y'  ,  I#1  / 


wn  I  J  \*i,  ;  ^  < 


Using  Eq .  (136),  the  result  is: 


2S+  A+  = 
in  m 


Y  fs+  +  s']  X 

/  -  l.  m  nj  mn 


Now,  using  Eq.  (126)  we  get: 


M2  £n< 


-2ikS  x 

z)f+(z)  e  n  dz 

m 


„  -2rkS  x 
2  m  .+ 

e  A 

m 


-ikfs++s”l x 

r+  r"  r+(z) £  (z)  e  m  n' 
m  n  m  n 

*  n 


r  T  X  e 
m  n  mn 


-ik.[s++S_  1 : 
l  m  nj 


and  it  follows  that: 


t*  E 


m  n=l 


lit  -  r  x  e 

c+  i.  mn 

1  Sm 


i.kfs+-s'|: 

(  m  nj 


Cons ider : 


f"*  i  dz 
m  n 


(141) 


(142) 


(143) 


where  A^  was  previously  defined  [see  Eq.  (47)].  Now  we  write: 


(144) 


(145) 


(146) 


(147) 


i 
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From  Eq .  (46)  we  now  have: 


f+  -$-=■  f  -  f  f+  +  k2  [fc'l2  -  fc+]2l  f+  f"  =  0  . 
m  az2  n  n  ^2  m  LI  nj  {  mj  J  m  n 


(148) 


However,  we  have  solved  equations  in  this  form  before,  so  it  follows: 
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f+  f"  dz 
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(149) 


Thus  the  connecting  formula  across  a  sharp  boundary  at  x  is: 
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2ik 
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S+  A+ 
m  m 


(150) 


RESULTS  FOR  TE. 

Results  for  the  TM  case  are  presented  in  Sec.  4,  in  this  subsec¬ 
tion,  we  discuss  some  results  for  the  TE  case.  In  the  TE  case,  we  are 
allowed  to  make  some  approximations  that  give  us  analytic  solutions 
for  some  important  quantities.  Unfortunately,  we  cannot  find  such 
solutions  in  the  TM  case  and  so  must  proceed  numerically.  First,  we 
consider  the  derivative  of  arid  C  with  respect  to  x,  recalling  that 


2 

qi 


1  +  C 


(151) 


It  follows  that 


qiqi.  *  CC' 


(152) 


Thus  , 
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v^v/AV'vy/v)! 


(153) 


Z!  =  -  -4  q!  -  -Z?  CC'  . 
1  2  ni  l 


To  find  C' ,  we  take  the  derivative  of  the  modal  equation  [Eq.  (95) 
and  find  that 


CC’\  ikh  +  Z 


if  K121 

i|  1  -  f 

-  / 


-  -  z'/z2 

c  6  g 


(154) 


Comparing  this  to  Eq .  (115)  shows: 


_i _ L__  z, 

2ik  A  CZ  g 
*  g 


(155) 


In  Eq.  (154),  the  (qg/C)^  dominates,  and  since  |  Zp  |  «  |kh|,  then 


£  C'  ikh  «  Z' 
C  g 


(156) 


We  integrate  this  directly  to  get 


C(x)  =  CQ  e 


ra  [Vx)-Zs,0,l 


(157) 


Equation  (157)  has  been  compared  to  the  numerical  solution  of  the 
modal  equation  [  Eq .  (95)]  at  5,  2.5,  and  50  kHz,  and  is  valid  to  less 
than  0.1  percent,  over  a  range  in  a g  from  4  S/m  down  to  10-;>  S/m. 

We  can  use  the  above  to  study  the  coupling  of  TE  modes.  From 
Eq .  (60)  it  is  clear  than  k^/K^  (where  m  *  n)  is  the  relative  amount 
of  mode  m  coupled  int o  mode  n  by  the  conductivity  change  from  x  to 
x  +  dx.  In  this  section  we  show  that  this  is  always  very  small,  and 
hence,  that  TE  mode  coupling  is  very  small,  as  expected. 

If  we  use  the  same  reasoning  as  above,  Eq .  (115)  can  be  ap¬ 
proximated  as 
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2 


<  158) 


then  Eq.  (122)  becomes,  for  the  TE  case  (recall  m  r*  n) 
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nm 
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kh 


The  diagonal  term  is 
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rm 
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So  we  have 


nn 


q'  C' 
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,q  C 
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Thus,  using  Eq .  (156)  we  get 
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(159) 
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This  approximation  appears  to  be  very  accurate  (greater  than  three 
decimal  places  at  20  kHz),  when  compared  to  the  complete  equation. 
Dividing  Eq .  (159)  by  Eq .  (162)  yields; 
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«  1  .  (163) 
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Thus,  the  of f-dlagonal  terms  are  very  small  and  we  do  not  expect  a 
large  amount  of  mode  coupling  for  the  TE  modes.  In  fact,  using 
Eq .  (157)  it  is  fairly  easy  to  show  that  Eq .  (70)  reduces  to 


Snm(x)  = 


2  n  1  S 

ikh  m2  -  n2  JS  (0)  S  (0)  q2 
m  n  g 
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Using  Eq .  (86),  we  find  that  1  4  Qn(x)  ~  1.  The  total  field, 


y*.  *>  «  £  r 

n-1 


>  Sn(x)  V  VX> 


(0)  -ik  J 

(x)  en(z’  X)  6 
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We  know,  however,  that 
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then,  is 

dx ' 
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(169) 
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so  it  follows  that 


VA  (°) 

e“(z’  i} 


V°>  Cn(i)  qR(i) 
V°)  qg(i)  Cn(^' 


sin  Ckh 


(170) 


en(z,  0) 


We  see,  therefore,  that  the  field  has  not  changed  significantly  across 
the  boundary. 

REFLECTIONS . 

In  this  subsection,  we  consider  the  reflection  term  bn  in 
Eq.  (58)  and  (59).  We  do  this  by  defining  a  reflection  coefficient  R 
such  that 
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R  a 
n  n 


(171) 


With  this  definition,  Eqs .  (58)  and  (59)  can  be  combined  into 
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where 
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(173) 


S  -  S 

-  n _ m 

®nm  S  +  S  ®nm 
n  m 


The  terms  in  the  summation  account  for  mode  coupling  and  are  clearly 
second  order.  We  have  shown  that,  at  least  for  TE,  they  can  be 
ignored.  Thus,  if  we  ignore  mode  coupling,  we  get: 


R(x)  = 


,,  -2ik  S  dx" 


(174) 


[This  form  gives  R(f)  =  0,  as  required,  since  there  can  be  no 
reflected  wave  past  x  *»  ^.]  To  determine  the  magnitude  of  R,  notice 
that  the  exponential  term  in  Eq.  (174)  is  less  than  one  due  to  the 
imaginary  part  of  S  (the  attenuation).  Since  we  do  not  expect  the 
attenuation  to  significantly  affect  the  wave  in  the  short  distance  of 
the  boundary,  we  write: 


R(x)  < 


dx'  , 


(175) 


which  we  solve  directly  as: 


R(x)  <  |  In  [S(i)/S(x)J 


(176) 


The  maximum  reflections  occur  at  x  *  0,  so, 


R(x.)  <  R(0)  <  2  fn  [S(f)/S(0) 


(177) 


Using  the  approximation  ior  the  TE  rnode--Lq.  (157)--we  can  show  (after 
some  manipulation) ,  that  if  the  change  in  Zg  across  the  boundary  is 
AZg ,  then 
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(178) 


R  “  TiT  AZ 
2kh  g 


Thus,  the  condition  for  negation  of  reflection  is: 


1  J_ 

2  kh 


AZ 

g 

so 

«  1 


(179) 


That  condition  is  met  in  the  TE  case  by  all  reasonable  boundaries. 

For  the  TM  case,  we  can  calculate  the  first-order  term  directly. 
The  reflection  term  |P.|  is  the  largest  for  the  Brewster's  mode  where 
C  ~  Zg.  For  a  frequency  of  30  kHz,  |R|  ~  0.02  and  changes  little  over 
the  VLF  range . 
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